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k*" ' Abstract 



We address the problem of parameter estimation for diffusion driven stochastic volatil- 
ity models through Markov chain Monte Carlo (MCMC). To avoid degeneracy issues we 
introduce an innovative reparametrisation defined through transformations that operate 
on the time scale of the diffusion. A novel MCMC scheme which overcomes the inherent 
difficulties of time change transformations is also presented. The algorithm is fast to 
, implement and applies to models with stochastic volatility. The methodology is tested 

■ through simulation based experiments and illustrated on data consisting of US treasury 

bill rates. 
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1 Introduction 

Diffusion processes provide natural models for continuous time phenomena. They are used 
extensively in diverse areas such as finance, biology and physics. A diffusion process is defined 
through a stochastic differential equation (SDE) 

dX t = n(t, X t , 9)dt + a(t, X t ,6)dW t , < t < T, (1) 
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where W is standard Brownian motion. The drift and volatility <r(.) reflect the instan- 
taneous mean and standard deviation respectively. In this paper we assume the existence 
of a unique weak solution to ([I]), which translates into some regulari ty conditions (locally 



Lipsc 



htz with a linear growth bound) on /x(.) and <r(.); see chapter 5 of 
19941 ) for more details. 



Rogers and Williams 



The task of inference for diffusion processes is 



particula r 



y cha llenging and has received 



Sorensenl (|2004l ) for an extensive review. 



remarkable attention in the recent literature; see 
The main difficulty is inherent in the nature of diffusions which are infinite dimensional ob- 
jects. However, only a finite number of points may be observed and the marginal likelihood 
of these observations is generally unavailable in closed form. This has sti mulated the devel- 



opme nt of various non-likeli hood approaches which use indirect inference (jGourieroux et al 



993T). estimating func t ions (IBibbv and Sorensen 



(Gallant and Tauchen 



1996); see also 



1995}, or the efficient method of moments 



Gallant and Lonel (119971 ). 



Most likelihood based methods approach the likelihood function through the transition 
density of ([I]). Denote the observations by Y k , k = 0, . . . , n, and with t k their corresponding 
times. If the dimension of equals that of X (for each k) we can use the Markov property 
to write the likelihood, given the initial point Yq, as: 



c(y, e\Y ) = n p k (Y k \Y k -r,e, A), A = t k - t k . 



k=l 



The transition densities p k (.) are not a vailable in c 



available. ' 
based, see 



hey may be ana l ytical, see 



Pedersen 



Ai't-Sahalia 



osed form but severa 



(2002) 



(j 19951 ). iDurham and Gallantl (|2002l ). They usually approximate the 



(2) 



app roximations are 



Ai't-Sahalial (120051 ). or simulation 



likelihood in a way so that t he discretisati o n err or can become arbitrarily small, although 
the methodology developed in lBeskos et al.1 (120061 ) succeeds exact inference in the sense that 



it allows only for Monte Carlo error. A potential downside of these methods may be their 
dependence on the Markov property. In many interesting multidimensional diffusion models 
the observation regime is different and some of their components are not observed at all. 



to model financial tim e series such as equity prices 



Stein and Steinl . 



1991 



ic volatility mo 


(Hcstorj 




1993 


: 



) , or interest rates (lAndersen and Lund 



1998; 



Hull and White , 



Durham! . 



2002; 



1987; 



Gallant and Tauchen 



199a). A stochastic volatility model is usually represented by a 2-dimensional diffusion 
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(3) 



dX t \ _ I n x {X t ,a t ,0) \ a x (a t ,9) \ / dB t 

da t J \ fi a (a t ,9) J \ cr a (a t ,9) J \ dW t 
where X denotes the observed equity (stock) log-price or the short term interest rate with 
volatility cr x (.), which is a function of a latent diffusion a. For the diffusion in ([3]), the Markov 
property may no longer hold; the distribution of a future stock price depends (besides the 
current price) on the current volatility which in turn depends on the entire price history. 
Stochastic volatility models are used 

An alternative approach to the problem adopts Bayesian inference utilizing Markov chain 
Monte Carlo (MCMC) methods. Adhering to the Bayesian framework, a prior p(9) is first 
assigned on the parameter vector 9. Then, g iven the observations Y, the posterior p(9\Y) 



can be explored through data augmentation (jTanner and Wong 



19871 ). treating the unob- 



served paths of X (paths between observations) as missing data. The resulting algorithm 



alternates bet ween updatin g 9 and X. Initia l MCMC schemes fol 

Jones ( 20031 



introduced by 



owing this programme wer e 



Erakei (2001) and 



Elerian et al 



Elerian et al. 



(I200lh . 



Jonesl (|1999l ); see also 
However, as not ed in the simulation ba s ed exp eriment of 

theoretically by iRoberts and Stramerl ()200ll ) , any such algorithm's convergence properties 



(|200lh and established 



will degenerate as the numbe r of imputed points increase s. The problem may be overcome 



with the reparametrisation of 



Roberts and Stramerl (|200ll ). and this scheme may be applied 



in all one-dimensional and some multi-dimensional contexts. However this framework does 
not co ver general multidimensional diffusion models. 



Chib et al 



(120051 ^ and 



Kalogeropoulos 



(|2007l ) offer appropriate reparametrisations but only fo r a class of stochastic volatility mod- 



els. Alternative reparametrisations were introduced in 



also 



Golightly and Wilkinson! (J2006|) for a sequential approach. 



Golightly and Wilkinson! (120071 ): see 



In this paper we introduce a novel reparametrisation that, unlike previous MCMC ap- 
proaches, operates on the time scale of the observed diffusion rather than its path. This 
facilitates the construction of irreducible and efficient MCMC schemes, designed appropri- 
ately to accommodate the time change of the diffusion path. Our approach is general enough 
to cover almost every stochastic volatility model used in practice. The paper is organized as 
follows: Section [2] elaborates on the need for a transformation of the diffusion to avoid prob- 
lematic MCMC algorithms. In Section [3] we introduce time change transformations whereas 
Section U] provides the details for the corresponding non-trivial MCMC implementation. The 
proposed methodology of the paper is tested and illustrated through numerical experiments 



in section [5l and on US treasury bill rates in section [6j Finally, section [7] concludes and 
provides some relevant discussion. 



2 The necessity of reparametrisation 

A Bayesian data augmentation scheme bypasses a problematic sampling from the posterior 
tt(9\Y) by introducing a latent variable X that simplifies the likelihood C(Y; X , 9). It usually 
involves the following two steps: 

1. Simulate X conditional on Y and 9. 

2. Simulate 6 from the augmented conditional posterior which is 
proportional to C(Y; X, 9)ir(9) . 

It is not hard to adapt our problem to this setting. Y represents the observations of the 
price process X. The latent variables X introduced to simplify the likelihood evaluations are 
discrete skeletons of diffusion paths between observations or entirely unobserved diffusions. 
In other words, X is a fine partition of multidimensional diffusion with drift fix(t, Xt, 9) and 
diffusion matrix 

E x (t,X t ,9) = a(t,X t ,9) x a(t,Xt,e)', 

and the augmented dataset is Xis, i = 0, . . . , T/S, where 5 specifies the amount of augmen- 
tation. The likelihood can be approximated via the Euler scheme 

T/5 

C E {Y;X,9) = f[p(X i5 \X (i _ 1)s ), X iS \X {i _ 1)s ~ M (X {i _ 1)s + Spt x (.), SE X (.)) , 

i=l 

, 199.4 ). 



which is known to converge to the true likelihood C(Y; X, 9) for small 5 (IPedersenl . 

Another property of diffusion processes relates to the quadratic variation process. 

Specifically we know that 



lim Y] { x iS ~ X(i-\)s) ~ x (i-x)s) = / ^x(s, X s ,9)ds a.s. 
i=i J o 



T/S T f T 

m ^ (Xt/; — Xu—-\\/i) {Xif, — = / 

The solution of the equation above determines the diffusion matrix parameters. Hence, 
there exists perfect correlation between these parameters and X as 5 — > 0. This has dis- 
astrous implications for the mixing and convergence of t he MCMC chain as it trans lates 



into reducibility for 5^0. This issue was first noted by 

scalar diffusions and also confirmed by the simulation experiment of 
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Roberts and Stramer 



2001 



Elerian et al 



for 



(2001 



Nevertheless, it is not an MCMC specific problem. It turns out that the convergence of its 
deterministic analogue, EM algorithm, is problematic when the amount of information in 
the augmented data X strongly exceeds that of the observations. In our case X contains an 
infinite amount of information for 5^0. 

The problem may be resolved if we apply a transformation so that the algorithm based 



Roberts and Stramerl (120011 ) 



on the transformed diffusion is no longer reducible as 5 — > 0. 
provide appropriate diffusion transformations for scalar diffusions. In a multivariate con- 
text t his requires a transformation to a diffusion with unit volatility matrix; see for in- 



stance 



Kalogeropoulos et al. 



((2007|). 



Ai't-Sahalia 



(120051 ) terms such diffusions as reducible 



and proves the non-reducibility of stochastic volatility models that obey ([3]). The transfor- 
mations introduced in this paper follow a slightly different route and target the time scale of 
the diffusion. One of the appealing features of such a reparametrisation is the generalisation 
to stochastic volatility models. 



3 Time change transformations 

For ease of illustration we first provide the time change transformation and the relevant 
likelihood function for scalar diffusion models with constant volatility. Nevertheless, one 
of the main advantages of this technique is the applicability to general stochastic volatility 
models. 

3.1 Scalar diffusions 

Consider a diffusion X defined through the following SDE: 

dX t = n(t, X t , 6)dt + adW t x , < t < 1 a > 0. (4) 

Without loss of generality, we assume a pair of observations Xq = y$ and X\ = y\ . For more 
data, note that the same operations are possible for every pair of successive observations 
that are linked together through the Markov property. We introduce the latent 'missing' 
path of X for < t < 1, d enoted by X mis , so that X = (y , X mis , yi ). In the spirit of 



Roberts and Stramerl (|200ll ). the goal is to write the likelihood for 9, a with respect to 
a parameter-free dominating measure. Using Girsanov's theorem we can get the Radon- 
Nikodym derivative between the law of the diffusion X, denoted by P , and that of the 
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driftless diffusion M = adW x which represents Wiener measure and is denoted by W x . We 
can write 



1 ^ = G(t,XM=exp[J o a2 dX s - - ^ - 2 dsj. 

By factorizing W x = W x x Leb(y{) x f(yi,a 2 ), where y\ ~ Af(yo,a 2 ) and Leb(.) denotes 



Lebesgue measure, we obtain 

dF(X mis ,y , yi ] 



G(t,X,0,a)f( yi] a), 



d{W x x Leb{y)} 

where clearly the dominating measure depends on a, since it reflects a Brownian bridge with 
volatility a. 

Now consider the time change transformation which first introduces a new time scale 
r](t, sigma)) 

<q(t,o) = [ a 2 ds = a 2 t, (5) 
J o 

and then defines the new transformed diffusion U as 



X v -i-(t,a)i 0<t<a 2 , 
M v -i itl a), t>a 2 . 

The definition for t > a 2 is needed to ensure that tU is well defined for different values of 
a 2 > which is essential in the conte x t of a MCMC algorithm. Using standard time change 
properties, see for example lOksendall (|2000l ). the SDE for U is 



dUt 



±tH(t, U t , 6)dt + dW t u 0<t<a 2 , 
dWY, t > a 2 , 



where W u is another Brownian motion at the time scale of U. By using Girsanov's theorem 
again, the law of U, denoted by P, is given through its Radon-Nikodym derivative with 
respect to the law of the Brownian motion W u at the [/—time scale: 



dF 
dW u 



G(t,U,9,a) = expj^ 



fj,(s,U s ,9) 



dU« 



fi(s,U s ,9f 



ds 



exp 



° ti(s,u a ,e) 



dU, - - 



i r n(s,u s ,ef 



ds 



(6) 
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If we condition the Wiener measure on y at the new time scale, the likelihood can be written 
with respect to a Brownian bridge measure as 



dP(U,y ,yi) = G(t,U,9,a)f( yi ;a)d{wV x Leb(y)} . 

However, this Brownian bridge is conditioned on the event U a i = y\ and therefore contains 
the parameter a. For this reason we introduce a second transformation which applies to both 
the diffusion's time scale and its path. Define 



U? = (a 2 -t)Z t/{a 2 {cT 2_ t)] , 0<t<a 2 , (7) 
t . t 



U? = U t -(l- —)y - —yi, 0<t<a 2 . 
Note that this transformation is 1-1. Its inverse is given by 

Zt = l+ J - U *H/(l+*n), 0<t< +00. 

Applying Ito's formula and using time change properties we can also obtain the SDE of Z 
based on another driving Brownian motion W z operating at the Z— time: 



dZ t = — i — ^ dt + dWf, 0<t< +00, (8) 

1 + ta z 

where v{Zt, a) = lit- This operation essentially transforms to a diffusion that runs from to 
+00 preserving the unit volatility. We can re-attempt to write the likelihood using Girsanov 
theorem and condition the dominating measure on y\ to obtain , 

d¥(Z,y , yi ) 

- G{t,Z,9,a)f{y l ;a). (9) 



d{Wf x Leb(y)} 

Despite the fact that G(Z,9,a) contains integrals defined in (0, +00), it is always finite 
being an 1 — 1 transformation of the Radon-Nikodym derivative between P(U) and given 
by ©. Using the following lemma, we prove that is the law of the standard Brownian 
motion and hence the likelihood is written with respect to a dominating measure that does 
not depend on any parameters. 

Lemma 3.1 Let W be a standard Brownian motion in [0, +00). Consider the process defined 
forO<t<T 

B t = (T- t)W t/mT - t )} + (1 - ^)yo + ^yi, < t < T 
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Then B is a Brownian bridge from yo at time to y\ at time T. 



Proof: See ([Rogers and Williamsl . Il994l . IV. 40.1) for the case yo = 0, T = 1. The extension 



for general yo and T is trivial. 

Corollary 3.1 The measure Wy is standard Wiener measure. 

Proof: Note that reflects a Brownian bridge from yo a t time to y\ at time T and 
we obtained by using the transformation of Lemma 13.11 Since this transformation is 
1 — 1, U is a Brownian bridge (under the dominating measure) if and only if Z is standard 
Brownian motion. 

Note that is the probability law of the driftless version of the conditional diffusion Z, 
whereas the SDE in ([8]) corresponds to the unconditional version of Z itself. The conditional 
SDE of Z is generally not available but this does not create a problem. For the path updates 
we may use the fact that 

P ,(Z\yo, yi) = G(t, Z, 9, <r) ^' ,<T \ oc G(t, Z, 9, a), (10) 



dwf li/u '* i; V ' ' ' 7 P (yi;a) 
where F y is the law of the conditional version of Z and f p (.) is the density of y\ under P. 
Both ¥ y and f p (.) are generally unknown but G(.) and /(.), which appear in ([9]) and (I10p . 
are available. 



3.2 Stochastic volatility models 

Consider the general class of stochastic volatility models with SDE given by ([3]). Without 
loss of generality we may assume a pair of observations (Xq = yo, X\ = y{) due to the 
Markov property of the 2-dimensional diffusion (X,a). The likelihood can then be divided 
into two parts: The first contains the marginal likelihood of the diffusion a and the remaining 
part corresponds to the diffusion X conditioned on the path of a 



F (X,a) =Fg(a)Fg(X\a). 
Denote the marginal likelihood for a by C a (a, 9). To ov ercome redu c tibilit y iss ues arising from 



the pa ths of a one may use the reparametrisations of 
(|2007l ). The relevant transformations of the latter are 



Chib et al 



([2005|) or 



Kalogeropoulos 



fa = h(a t ,9), 



dh(ctt, 9) 
da t 



{a a (a t ,9)} 1 , 



7t = fit ~ Po, Pt = v(lt), 
and the marginal likelihood for the transformed latent diffusion 7 becomes 



By letting a.% = gj = h 1 (r)('yt), 9), the SDE of X conditional on 7 becomes: 

dX t = fi x {X t ,9l 6)dt + a x (gl 9)dB t , 0<t<l. 

Given the paths of the diffusion a, the volatility function cr x (g^,9) may be viewed as a 
deterministic function of time. The situation is similar to that of the previous section. We 
can introduce a new time scale 

#,7,0) = f 4(gl9)ds, 
Jo 

and define U with the new time scale as before (M is a Brownian motion on the [/—time 
scale) 



Xri-lff), < t < T, 

U t ={ v (} - - (12) 



M v - Ht) , t>T. 



The SDE for U now becomes 



dUt = {^p^n^\ dt + dw v^ 0<t<T . 

[ °"x(7^- 1 (t, 7 ,0) 1 c') J 

We obtain the Radon Nikodym derivative between the distribution of U with respect to that 
of the Brownian motion W u , 

d¥ 

G(U,j,9), 



dW u 

and introduce as before. The density of y\ under W u , denoted by f(y, 7, 9), is just 

f( yi ;y,9) = N(y ,T). 

The dominating measure reflects a Br ownian motion conditioned to equal y at a 
parameter depended time T = rj(tk+i, 7, 9). To remove this dependency we introduce a 
second time change 

9 



U? = (T — t)Z t/{T(T _ t)} , < t < T, (13) 
U° = U t -(l- ^)y - ^yi, < t < T. 

Therefore, the SDE for Z is now given by 

dz * = TT^ m L + v(Zt)\dt + dW?, 0<t<^, 

1 + tT ( 0£(7fc(t,7,«)>0) 

where k(t,^,9) denotes the initial time scale of X and v{Zt) = Ut- 

Conditional on 7, the likelihood can be written in a similar manner as in ([9]): 

(Z\y ,y ul )=G(Z, 1 ,e)f(y 1 ; 1 ,e) (14) 



d{Wf x Leb{y)} 



It is not hard to see that reflects a standard Wiener measure and therefore the dominating 
measure is independent of parameters. To obtain the full likelihood we need to multiply the 
two parts given by (jlip and (|14p . 

3.3 Incorporating leverage effect 

In the previous section we made the assumption that the increments of X and 7 are inde- 
pendent, in other words we assumed no leverage effect. This assumption can be relaxed in 
the following way: In the presence of a leverage effect p, the SDE of X conditional on 7 can 
be written as 



dX t = p x (X t ,gl 9)dt + pa x (g], 0)dW t + y/l - p 2 a x (gl 6)dB t , 0<t<t k , 

where W is the driving Brownian motion of 7). Note that given 7, W can be regarded as 
a function of 7 and its parameters 8. Therefore, the term pa x {gj ,9)dWt can be viewed as 
a deterministic function of time, and it can be treated as part of the drift of Xt. However, 
this operation introduces additional problems as the assumptions ensuring a weakly unique 
solution to the SDE of X are violated. To avoid this issue we introduce the infinitesimal 
transformation 

X t = H(H t ,p, 1 ,6)=H t + I pa x (g], 8)dW s , 

Jo 
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which leads us to the following SDE for H: 



dH t = fi x {H(X t ,p, 7, 0), gl,0} dt + ^l- p 2 a x (gl 9)dB t , 0<t< t k . 

We can now proceed as before, defining U and Z based on the SDE of H in a similar manner 
as in (fT2"j) and (fT3|) respectively. 

3.4 State dependent volatility 

Consider the family of state dependent stochastic volatility models where conditional on 7, 
the SDE of X may be written as: 

dX t = p x {X t ,gl 6)dt + ax{gl d)a 2 {X t , 0)dB t , 0<t< t k . 



his c l ass contains amo n g others , the m odels of 



Andersen and Lund 



(|1998h 



Gallant and Tauchen 



(|1998l ) , iDurhaml (|2002l ) , lErakerl (|200ll ) . In order to apply the time change transformations 
of section \3/2\ we should first transform X to Xt, through Xt = h(Xt,9), so that it takes 
the form of (El). Such a tran sformation, which may be viewed as the first transformation in 
Roberts and Strainer! (|200ll ). should satisfy the following differential equation 



dh(X t 



1 



dX t a 2 (X t ,9Y 

The time change transformations for U and Z may then be defined on the basis of X that 
will now have volatility o\(g1,9). The transformation h(.) also applies to the observations 
(2)0 1 Vi) which are now functions of parameters. This would translated in a parameter 
dependent likelihood dominating measure, had it n ot been for the second step i n (j 12 j) which 



Roberts and Stramer (2001 



). Note that 



in this case acts like the second transformation in 
the parameters of G2{Xt,6) enter the reparametrised likelihood in two ways: first through 
the f(y; 7, 6) which now should include the relevant Jacobian term, and second through the 
drift of Z which is centered at based on the transformed observations. 



3.5 Multivariate stochastic volatility models 



We may use the techniques of section 13.31 to define time change transformations for multidi- 
mensional diffusions. Consider a d— dimensional version of the SDE in flU) where a now is a 
2x2 matrix ([cr]ij = <Tij). As noted in 



Kalogeropoulos et al 



(|2007l ). the mapping between 



a and the volatility matrix aa T should be 1-1 in order to ensure identifiability of the a 
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parameters. A way to achieve this, is by allowing a to be the lower triangular matrix that 
produces the Cholesky decomposition of aa T . For d = 2, the SDE of such a diffusion is given 
by 

dx\ l} = n(x{ 1] , x/ 2} ,9)dt + andB t , 
dx} 2} = p(x} 1} ,x} 2} ,6)dt + a 21 dB t + u 22 dW t . 

The time change transformations for X^ will be exactly as in section [370 For X^ 
note that given X^ the term a 2 \dBt is now a deterministic function of time and may be 
treated as part of the drift. Thus, we may proceed following the route of the section [3~3l 

Similar transformations can be applied for diffusions that have, or may be transformed 
to have, volatility functions independent of their paths. For example we may assume two 
correlated price processes with correlation p x : 

[a] 21 =p x ai 2 \gl9), 

[a] 22 = ^7 x a x 2} (9l0). 

We may proceed in a similar manner for multivariate stochastic volatility models of general 
dimension d. 

4 MCMC implementation 

The construction of an appropriate data augmentation algorithm involves several issues. The 
time change transformations introduce three interesting features to the MCMC algorithm 
which we address separately: the presence of three time scales; the need to update diffusion 
paths that run from to +oo; and the fact that time scales depend on parameters. For 
ease of illustration we will assume the simple case of a univariate diffusion with constant 
volatility and a pair observations (Xq = yo and X\ = y\). Extensions and generalisations of 
the algorithm for stochastic volatility models are noted where appropriate. 

4.1 Three time scales 

We introduce m intermediate points of X at equidistant times between and 1, to give 
X = {Aj/( m +i), i = 0, 1, . . . ,m + 1}. In addition, we make the assumption that m is 
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large enough for accurate likelihood approximations and any error induced by the time 
discretisation is negligible for the purposes of our analysis. 

Given a value of the time scale parameter a, we can get the U— time points by applying 
© to each one of the existing points X, so that 

U a 2i/(m+l) = X i/{m+l), i = 0, 1, • • • ,m + 1. 

Note that it is only the times that change, the values of the diffusion remain intact. In a 
stochastic volatility model we would use the quantities 

i+l 

a 2 x {.)ds 

m+l 

for each pair of consecutive imputed points. 

The points of Z are multiplied by a time factor which corrects the deviations from unit 
volatility. The Z— time points may be obtained by 

7 a 2 i/(m + 1) „ _ 

t z = L± ! I 7=01 m 

% 111 2 • If i 1 \\ > ' > ' ' ' ' 

Clearly this does not apply to the last point which occurs at time +oo. Therefore, the paths 
of X, or U, are more convenient and may be used for likelihood evaluations exploiting the 
fact that the relevant transformations are 1-1. However, the component of the relevant Gibbs 
sampling scheme is the diffusion Z. 

Figure [1] shows a graphical representation of X, U and Z plotted against their corre- 
sponding time scales for a = y/2 and m = 7. Although X and U have the same values, their 
volatilities are v2 and 1 respectively. The ending point of Z does not appear on the graph 
as it occurs at time +oo. 

[Figure 1 about here.] 
4.2 Updating the paths of Z 

The paths of Z may be updated using an independence sampler with the reference measure 
as a proposal. Here W z reflects a Brownian motion at the Z— time which is fixed given 
the current values of the time-scale parameter (s). An appropriate algorithm is given by the 
following steps. 

• Step 1: Propose a Brownian motion on the Z— time, say Z* . The value at 
the endpoint (time +oo) is not needed. 
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• Step 2: Transform back to U* , using ([7]). 

• Step 3: Accept with probability: min |l, } ■ 

4.3 Updating time scale parameters 

The updates of parameters that define the time scale, such as a, are of particular interest. 
In almost all cases, their conditional posterior density is not available in closed form, and 
Metropolis steps are inevitable. However, the proposed values of these parameters will imply 
different Z— time scales. In other words, for each potential proposed value for a there 
exists a different set of Z— points needed for accurate approximations of the likelihood the 
Metropolis accept-reject probabilities. In theory, this would pose no issues had we been able 
to store an infinitely thin partition of Z, b ut of course this is not possible . 



We use retrospective sampling ideas; see 



Papaspiliopoulos and Roberts! (|2005l ) and 



Beskos and Roberts 



(2()()5j for applications in different contexts. Under the assumption of a sufficiently fine par- 
tition of Z, all the non-recorded intermediate points contribute nothing to the likelihood 
and they are irrelevant in that respect; the set of recorded points is sufficient for likelihood 
approximation purposes. Alternatively, we may argue that their distribution is given by the 
likelihood reference measure which reflects a Brownian motion. Thus, they can be drawn af- 
ter the proposal of the candidate value of the time scale parameter. To ensure compatibility 
with the recorded partition of Z, it suffices to condition on their neighboring points. This is 
easily done using standard Brownian bridge properties: Suppose that we want to simulate 
the value of Z at time % which fall between the recorded values at times t a and t c , so that 
ta < if) < t c . Denote by Zt a and Zt c the corresponding Z values. Under the assumption that 
Z is distributed according to between t a and t c we have that 

z, t \ ~n { {tt ~ u)z '; + { !< ~ th)Zt - , ttb - l ; ){t ; - tb> ). us) 

The situation is pictured in Figure [21 where the black bullets represent the stored points and 
the triangles the new points required for a proposed value of a. The latter should be drawn 
retrospectively given the former via (|15p . 

[Figure 2 about here.] 

A suitable algorithm for the a— updates may be summarized through the following steps: 



• Step 1: Propose a candidate value for a, say a*. 
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• Step 2: Repeat for each pair of successive points: 

— Use d5j) and ([7]) to get the new times associated with it. 

— Draw the values of Z at the new times using (1151) . 

— Transform back to U* , using (JTD . 

Form the entire path U* by appropriately joining its bits. 

• Step 3: Accept with probability: min 1 1, G ^jj* 'q J^f^y^ ^ } • 

Note that in a stochastic volatility model the paths of the transformed diffusion jt are 
associated with the time scale of the Z. Therefore a similar algorithm may be used for their 
updates. 

5 Simulations 

As discussed in section [2j appropriate reparametrisations are necessary to avoid issues re- 
garding the mixing and convergence of the MCMC algorithm. In fact, the chain becomes 
of augmentation incr eases. This is also verified by the numerical ex- 



reducible as the leve 

amples performed in iKalogeropoulosI (|2007l ) even in very simple stochastic volatility models. 
In this section we perform a simulation based experiment to check the immunity of MCMC 
schemes to increasing levels of augmentation, as well as the ability of our estimation proce- 
dure to retrieve the correct values of the diffusion parameters despite the fact that the series 
is partially observed at only a finite number of points. We simulate data from the following 
stochastic volatility model 



dX t = K x (fi x - X t )dt + pexp(a t /2)dW t + y/l - p 2 exp(a t /2)dB t , 
dat = K a (fji a — atjdt + adWt, 



where B and W are independent Brownian motions, and p reflects the correlation between the 
increments of X and a, also term as leverage effect. A high frequency Euler approximating 
scheme with a step of 0.001 was used for the simulation of the diffusion paths. Specifically, 
500, 001 points were drawn and one value of X for every 1000 was recorded, thus forming a 
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dataset of 501 observations of X at < t < 500. The parameter values were set to p = —0.5, 
a = 0.4, k x = 0.2, p x = 0.1, K a = 0.3 and p a = —0.2 

The transformations required to construct an irreducible data augmentation scheme are 
listed below. First we transform a to 7 through 



ott — an 

It = — -, < t < 500, 

a 



Given 7, and for each pair of consecutive observation times t^-i and (k = 1,2,..., 500) 
on X, we transform as follows: First, we remove the term introduced from the leverage effect 

H t = X t - / pexp {v(~f s , a, a )/2} dW s , t k _ x <t<t k , 
Jtk-i 

and consequently we set 

V(t) = (! - p) 2 exp{v(j s ,a,a }ds. 
Jtk-i 

Then, U and Z may be defined again from [12] and [13] respectively, but based on H rather on 
X. The elements of the MCMC scheme are Z,*y, an and the parameters (k x , fi x , K a , fi a , p, a). 

We proceed by assigning flat priors to all the parameters, restricting k x , K a , a to be 
positive and p to be in (—1, 1). The number of imputed points was set to 30 and 50, the length 
of the overlapping blocks needed for the updates of 7 was 2, and the relevant acceptance 
rate 75% whereas the acceptance rate for X was 95%. Figure [3] shows autocorrelation plots 
for the 2-dimensional diffusion's (X, a)' volatility parameters p and a. There is no sign of 
any increase in the autocorrelation to raise suspicions against the irreducibility of the chain. 
Figure H] shows density plots for all parameters and both values of m. These plots indicate a 
sufficiently fine discretisation and a good agreement with true values of the parameters. The 
latter is also confirmed by Table [1] 

[Figure 3 about here.] 
[Figure 4 about here.] 
[Table 1 about here.] 
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6 Application: US treasury bill rates 



To illustrate the time change methodology we fit a stochastic volatility model to US treasury 
bill rates. The dataset consists of 1809 weekly observations (Wednesday) of the 3— month 
US Treasury bill rate from the 5th of January 1962 up to the 30th of August 1996. The data 
are plotted in Figure G3 

[Figure 5 about here.] 



Pre yious analyses of these data include 



Durham! ( 



2002) 



Andersen and Lundl (119981 ) 



Durham and Gallantl ( 



20021 ) 



Eraker 



2001 



Gallant and Tauchen 



, and 



Golightly and Wilkinson 



20061 ). Apart from some slight deviations the adopted stochastic volatility models consisted 



of the following SDE. 



dr t = (0 o -0 1 r t )dt + rtexp(a t /2)dBt, 
dat = k(h — ottjdt + adWt, 



(16) 



with independent Brownian motions B and W . In some cases the following equivalent model 
was used: 



dr t = (0q - 0\r t )dt + a r rf exp(a t /2)dB t , 
dat = —Kottdt + adWt- 



(17) 



We proceed with the model in ([16 ), as posterior draws of its p aram eters exhibit substantially 



l ess au tocorrelation. In line with 
(120061 ). we also set ip = 1 



Gallant and Tauchenl (119981 ) and 



Eraker ( 



2001 



Durham ( 



2002 ) and 



Golightlv and Wilkinson 



Durham and Gallant 



] ( 20021 



assume general 'elasticity of variance' ip but their estimates do not indicate a significant 
deviation from 1. By setting Xt = log(rt), the volatility of Xt becomes exp(aj/2). Therefore 
the [/—time for two consecutive observation times tk-i and i& is defined as 



v(t) 



exp(at)ds, 



and U and Z are given by (|12p and (|13p respectively . We also transform a to 7 as before: 



It 



at - a 



a 



17 



at = "(it, oco) = «o + <?Jt- 

We applied MCMC algorithms based on Z and 7 to sample from the posterior of the 
parameters 9q, $i> K > A* an d cr. The time was measured in years setting the distance between 
successive Wednesdays to 5/252. Non-informative priors were assigned to all the parameters, 
restricting k and a to be positive to ensure identifiability and eliminate the possibility of 
explosion. The algorithm was run for 50, 000 iterations and for m equal to 10 and 20. To 
optimize the efficiency of the chain we set the length of the overlapping blocks of 7 to 10 
which produced an acceptance rate of 51.9%. The corresponding acceptance rate for Z was 
98.6% . 

The kernel density plots of the posterior parameters and likelihood (Figure [6]) indicate 
that a discretisation from an m of 10 or 20 provide reasonable approximations. The cor- 
responding autocorrelation plots of Figure [7] do not show increasing autocorrelation in m, 
a feature that would reveal reducibility issues. Finally, summaries of the posterior draws 
for all the parameters are provided in Table [2J The parameters k, (i and a are different 
from verifying the existence of stochastic volatility. On the other hand, there is no evi- 
dence to support the existence of mean r eversion o n the r ate, as 9g and 9-\ are not far from 



0. The results are in line wit 



1 those of 



Golightlv and Wilkinson! pOOfl ). 



Durham! ([2002), Durham and Gallant! ((2002) and 



7 Discussion 



[Figure 6 about here. 



[Figure 7 about here. 



[Table 2 about here.] 



Data augmentation MCMC schemes constitute a very useful tool for likelihood-based infer- 
ence on diffusion models. They may not have the appeal ing properties of complete elimi- 



nation of the time discretisation error 
likelihood expressions of 



Bcs 



cos et al 



20061 ), or the closed form approximate 



Ai't-Sahalial ([2002]), but nevertheless they give a satisfactory and 
very general solution to the problem. However data augmentation schemes require careful 
construction to avoid the degeneracy issues described at the beginning of this paper. 

Here, we introduce an innovative transformation which operates by altering the time axis 
of the diffusion. To accommodate the special features of time change transformations we also 
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introduce a novel efficient MCMC scheme which mixes rapidly and is not provibitively com- 
putationally expensive. Our method is also easy to implement and introduces no additional 
approximation error other than that included in methodologies based on a discretisation of 
the diffusion path. Moreover it is general enough to include general stochastic volatility 
models. 

Further work will consider problems with state-dependent volatility and models which 
involve jump diffusions, to which the methodology introduced here can be easily applied. 
Fundamental to our approach here has been the introduction of a non-centered parameter- 
isation to decouple dependence inherent in the model between missing data and volatility 
parameters. However non-centered constructions are not unique, a s illustrated by the choice 



i n the diffusion context betw e en th e state rescaling approaches of 



Golightlv and Wilkinson 



20071 ): 



Roberts and Stramerl (|200ll ) and the time-stretching strategy adopted here. Clearly, 
further work is required to investigate the relative merits of these approaches in different 
situations. 
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Figure 1: Plots of a sample path for X, U and Z against their corresponding times for 
a = y/2 and m = 7. Z equals at time +oo. 
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Figure 2: Updates of time scale parameters: For every proposed value of them, new points 
are required and should obtained conditional on the stored points. 
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Figure 3: Autocorrelation plots for the posterior draws of p and a for different numbers of 
imputed points m = 30, 50. Simulation example of Chapter 3. 
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Figure 4: Kernel densities of the posterior draws of all the parameters for different numbers 
of imputed points m = 30, 50. Simulation example of Chapter 3. 
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Figure 5: Weekly 3— month US Treasury bill rate from the 5th of January 1962 up to 
30th of August 1996. 
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Figure 6: Kernel densities of the posterior draws of all the parameters and the log-likelihood 
for different values of imputed points m = 10, 20. Example on Weekly 3— month US Treasury 
bill rates. 
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Figure 7: Autocorrelation plots for the posterior draws of the model parameters for different 
numbers of imputed points m = 10, 20 for the analysis of Weekly 3— month US Treasury bill 
rates. 
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Parameter 


True value 


Post, mean 


Post. SD 


Post 2.5% 


Post median 


Post 97.5% 




0.2 


0.244 


0.038 


0.173 


0.243 


0.321 




0.1 


0.313 


0.174 


-0.046 


0.317 


0.641 


K a 


0.3 


0.304 


0.148 


0.110 


0.277 


0.672 




-0.2 


-0.268 


0.107 


-0.484 


-0.267 


-0.059 


a 


0.4 


0.406 


0.130 


0.202 


0.390 


0.705 


P 


-0.5 


0.477 


0.138 


-0.657 


-0.491 


-0.066 



Table 1: Summaries of the posterior draws for the simulation example of Chapter 3 for 
m = 50. 
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Parameter 


Post, mean 


Post. SD 


Post 2.5% 


Post median 


Post 97.5% 




0.130 


0.238 


-0.347 


0.132 


0.589 


01 


0.013 


0.057 


-0.096 


0.013 


0.125 


K 


2.403 


0.620 


1.319 


2.360 


3.745 




-3.966 


0.211 


-4.384 


-3.964 


-3.547 


a 


2.764 


0.311 


2.199 


2.750 


3.420 



Table 2: Summaries of the posterior draws for the stochastic volatility model of Weekly 
3— month US Treasury bill rates. 
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